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SUMMARY 


This  paper  presents  an  extension  of  a  structural  finite  element  mesh 
improvement  technique  to  heat  conduction  analysis.  The  mesh  improvement  con¬ 
cept  was  originally  presented  by  Prager  in  studying  tapered,  axially  loaded 
bars.  It  was  further  developed  by  Kittur,  et  al . ,  who  showed  that  an  improved 
mesh  can  be  obtained  by  minimizing  the  trace  of  the  stiffness  matrix.  In  this 
paper  these  procedures  are  extended  and  applied  to  the  analysis  of  heat  conduc¬ 
tion  in  an  infinitely  long  hollow  circular  cylinder. 
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NOMENCLATURE 

shape  function  for  element  e 
number  of  elements 
temperature  gradient,  °C/mm 
radius,  mm 

inner  radius  of  the  cylinder,  mm 
outer  radius  of  the  cylinder,  mm 
inner  radius  of  element  e,  mm 

inner  radius  of  element  e+1  or  outer  radius  of  element  e,  mm 
temperature  at  the  inner  radius,  °C 
temperature  at  the  outer  radius,  °C 
stiffness  matrix  or  element  e 
thermal  conductivity,  W/mm  °C 
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thermal  conductivity  of  the  material  of  element  e,  W/mm  °C 
trace  of  the  global  stiffness  matrix 
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INTRODUCTION 

The  finite  element  analyst  is  continually  confronted  with  the  decision 
of  selecting  a  "good"  mesh  for  a  given  problem.  Unfortunately  there  seems  to 
be  no  simple  mesh  optimization  method  which  is  generally  applicable.  Recently 
Kittur,  et  al.,  (ref.  1)  have  proposed  a  method  which  provides  an  improved  mesh 
over  a  uniform  mesh  for  a  broad  class  of  structural  problems.  The  method  is 
based  upon  the  simple  procedure  of  node  relocation  to  minimize  the  trace  of  the 
stiffness  matrix.  Figure  1  presents  a  flowchart  illustrating  the  procedure. 

An  advantage  of  this  method  is  that  the  analyst  may  select  the  nodal 
positions  prior  to  solving  the  equilibrium  equations.  A  posteriori  methods 
such  as  those  of  Shephard  (ref.  2)  and  Carroll  (ref.  3)  can  then  be  used  to 
refine  the  mesh  further. 

Kittur,  et  al . ,  have  shown  that  the  minimum  trace  method  produces 
excellent  results  for  the  problem  studied  by  Prager  (ref.  4):  a  uniformly 
tapered,  axially  loaded  bar.  This  problem  is  governed  by  Poisson's  equation 
with  Dirichlet  boundary  conditions. 

In  this  paper  we  extend  and  apply  the  method  to  analyze  the  heat  con¬ 
duction  through  the  cross  section  of  an  infinitely  long  hollow  circular  cylin¬ 
der.  For  this  problem,  the  governing  equation  is  also  of  the  form  of  Poisson's 
equation.  In  this  case,  however,  since  no  heat  sources  are  considered  the 
governing  equation  reduces  to  Laplace's  equation.  Both  Dirichlet  and  mixed 
boundary  conditions  are  considered.  The  results  are  compared  with  analytical 
solutions  and  with  finite  element  solutions  obtained  with  a  uniform  mesh. 


CONFIGURATION  AND  PROBLEM  DEFINITION 


Consider  a  hollow  circular  cylinder  of  infinite  length  having  inner  and 
outer  radii:  r0  and  rn.  Let  the  thermal  conductivity  be  k,  and  the  temper¬ 
atures  at  the  inner  and  outer  radii  be  T0  and  Tn  respectively.  The  govern¬ 
ing  equation  for  the  temperature  distribution  along  a  radial  line  is: 


The  boundary  conditions  are: 


and 


T  =  T  at  r  =  r  , 
o  o 

T  =  T  at  r  =  r 
n  n 


The  solution  of  equation  (1)  subject  to  equation  (2)  is: 


T*Tn* 
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Next,  suppose  that  the  temperature  gradient  at  the  inner  surface  is 
specified  as:  q0.  The  boundary  conditions  are  then: 


and 


dT 
dr 

T  =  T. 


dr  =  qo 


n 


at  r  =  r 


at  r  =  r. 


In  this  case  the  solution  of  equation  (1)  is: 


T  =  T  -  r  q  1  n  — 
n  oMo  \ r 


FINITE  ELEMENT  FORMULATION  Af'D  MESH  OPTIMIZATION 


(4) 


(5) 


Figure  2  shows  the  finite  element  mode1.  Its  consists  of  a  series  of 
annular  elements.  For  element  (e)  let  the  inner  and  outer  radii  be  re  and 
re+].  The  entries  of  the  element  stiffness  matrix  a:*.-: 


dN1:  dN! 


r  dT  dT  dr 


(6) 


where  k  is  the  element  conductivity  constant  and  where  the  element  shape 
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functions  N]  and  N2  are: 
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By  carrying  out  the  indicated  operations  the  element  stiffness  matrix 
becomes: 
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where  is  defined  as: 
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Se  '  T,Ke\  r 


r  ,  +  r  ' 
e+1  e 


e+1  "  re/ 


(8) 


(9) 
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Hence  the  trace  x  of  the  global  stiffness  matrix  is: 


x  =  2 


where  n  is  the  number  of  elements. 

The  trace  may  be  minimized  with  respect  to  the  nodal  coordinates  by  set¬ 
ting  the  partial  derivative  of  x  with  respect  to  re  equal  to  zero  and  solv¬ 
ing  for  the  ratio  re+]/re: 


Since  the  derivative  of  a  sum  is  the  sum  of  the  derivatives,  and  assuming  that 
the  conductivity  k  is  uniform  thorugh  the  body  (i.e.:  k6_]  =  Ke  =  Ke+i>, 

(11)  becomes: 


L  re  *  re-l  re+l  *  re 

el  +  re  “  re-l  +  re+l  "  re  +  ' 


which  simplifies  to: 


r„  -  r„ 
>.  e  e 


re+l  "  re, 


Inverting  (13)  and  factoring: 


2  re 


e-1  r 


Rearranging: 


re  Yre+1 


re-l A  re 


re+l  1  {  re 


re  J  Vre-1 


Let  A  =  — —  and  B  -  - —  ,  then  (15)  becomes: 
re  re-l 


B(A  -  1)2  =  A(B  -  l)2 


AB(A  -  B)  -  (A  -  B)  =  0 
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thus: 


(A  -  BXAB  -  1)  =  0  (18) 

Equation  (18)  implies  that  either  (A  -  B)  =  0  or  (AB  -  1)  =  0.  If  (AB  -  1)  =  0, 

r_ . ,  r.  r.  , 


then  AB  =  1 ,  which  means 


re+1  re 
re  re-l 


which  is  impossible.  Therefore 


(A  -  B)  =  0,  which  implies  that  A  =  B.  This  leads  to  the  relatively  simple 
relation: 


where  y  is  a  constant.  The  optimal  nodal  positions  are  found  through 
repeated  use  of  equation  (19): 
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thus:  rn  =  r0yn  and  y  =  (rn/rg)  '  •  Changing  n  to  e,  we  have: 


re  =  r0Y  or: 


re  =  r0  V 


NUMERICAL  EXAMPLES 


Example  1 

To  illustrate  the  effectiveness  of  the  method  consider  an  annular  cylinder 
with  the  following  temperatures  specified  on  the  boundaries  and  made  from  a 
material  with  constant  thermal  conductivity  k  =  Ke  =  1.0. 

T  =  T0  =  100  °C  at  rn  =  20  mm 


T  =  Tn  =  0  °C  at  rn  =  50  mm 
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Consider  two  finite  element  models,  each  with  four  linear  elements  and 
five  nodes:  The  first  model  (designated  "Uniform  Mesh")  has  equally  spaced 
nodes.  The  second  model  (designated  "Improved  Mesh")  has  its  nodal  spacing 
governed  by  equation  (23).  The  objective  is  the  determination  of  the  internal 
temperature  distribution. 

The  solutions  of  the  finite  element  governing  equations  lead  to  the 
results  listed  in  table  I.  The  error  at  the  nodes  is  defined  as  the  difference 
between  the  theoretical  results  and  the  finite  element  results.  In  order  to 
compare  the  error  distribution  among  the  elements,  rms  errors  were  calculated 
using  50  uniformly  spaced  points  through  the  thickness  of  each  element.  Each 
element's  linear  shape  function  was  used  to  compute  the  temperature  at  any 
point  along  the  element  from  nodal  temperature  values.  Error  values  were 
calculated  at  each  point  from  the  difference  between  theoretical  and  finite 
element  temperature  values.  The  rms  errors  are  shown  beside  the  nodal  errors 
in  table  I.  The  standard  deviation  of  the  wi thin-element  rms  errors  is  shown 
at  the  bottom  of  table  I  with  the  overall  rms  error  for  the  model  which  was 
calculated  using  200  uniformly-spaced  points  through  the  model.  (These  are  not 
quite  the  same  points  used  for  computing  the  within-element  rms  errors.) 

The  theoretical  solution  for  this  example  (from  eq.  3)  is  shown  in 
figure  3.  The  nodal  values  for  both  models  are  superimposed  on  the  theoretical 
solution.  At  this  scale,  differences  between  nodal  and  FEM  values  are  much 
too  small  to  be  seen  in  the  figure.  The  errors  for  the  two  models  are  plotted 
on  an  expanded  scale  in  figure  4.  Since  linear  elements  were  used,  the  errors 
tend  to  increase  towards  the  center  of  each  element  (away  from  the  nodes). 

The  uniform  mesh  has  relatively  small  errors  at  the  nodes  and  greater 
errors  near  the  inner  surface  of  the  annular  cylinder  where  the  temperature 
gradient  is  greatest.  The  improved  mesh  is  found  to  have  zero  errors  at  the 
nodes  and  a  uniform  distribution  of  errors  between  elements.  The  overall  rms 
error  of  the  improved  mesh  is  23  percent  less  than  that  of  the  uniform  mesh. 


Example  2 


Next,  consider  the  same  cylinder  but  let  the  temperature  gradient  be  spec¬ 
ified  on  the  inner  boundary.  Specifically,  let  the  boundary  conditions  be: 


=  -5.4567833  °C/mm  at  rQ  =  20  mm 

T  =  T  =  0  °C  at  r  =  50  mm 

n  n 


(25) 


(The  boundary  conditions  were  chosen  to  yield  a  nearly  identical  theoreti¬ 
cal  temperature  distribution  as  in  the  first  example.) 

Table  II  shows  the  comparison  between  the  finite  element  solution  and  the 
theoretical  values  of  the  temperature.  Once  again  the  improved  mesh  shows  a 
more  uniform  distribution  of  rms  errors  between  elements  and  the  overall  rms 
error  is  lower  (by  15  percent)  than  that  of  the  uniform  mesh.  The  errors  for 
the  two  models  are  plotted  in  figure  5. 
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Example  3 


Since  the  temperature  gradient  is  specified  at  the  inner  boundary,  it  is 
also  of  interest  to  know  how  the  values  of  the  temperature  gradients  obtained 
using  the  two  finite  element  meshes  compare  with  each  other  and  with  the  theo¬ 
retical  values,  table  III  provides  such  a  comparison.  The  errors  as  plotted 
in  figure  6  are  much  greater  (relatively)  than  in  the  previous  examples,  espe¬ 
cially  at  the  nodes.  Note  that  the  largest  error  is  at  the  inside  radius  where 
the  derivative  boundary  condition  is  applied.  Generally  the  Neumann  boundary 
conditions  are  not  strongly  satisfied  in  the  finite  element  method  (ref.  5). 
Therefore,  in  general,  the  FE  method  is  not  well  suited  for  gradient-type  prob¬ 
lems.  The  improved  mesh  shows  a  slightly  more  uniform  distribution  of  rms 
errors  between  elements  and  the  overall  rms  error  is  also  slightly  lower  (by 
7  percent). 


DISCUSSION 

The  numerical  examples  show  that  the  "improved"  mesh  defined  by  equation 
(23)  produces  results  which  are  in  better  agreement  with  the  theoretical  values 
than  those  obtained  using  the  uniform  mesh.  Therefore,  the  mesh  defined  by 
equation  (12)  is  an  improvement  over  the  uniform  mesh  for  both  the  temperature- 
fixed  boui.dary  conditions  and  for  the  mixeu  boundary  conditions. 

In  Example  1,  the  values  of  the  stiffness  matrix  traces  of  the  uniform  and 
improved  meshes  are  74.6667  and  70.1520,  respectively.  This  indicates  that  for 
this  example  the  trace  is  only  moderately  sensitive  to  the  nodal  locations. 

More  dramatic  differences  in  the  results  will  occur  in  problems  where  the  trace 
is  more  sensitive  to  changes  in  the  nodal  positions.  This  will  occur  in  prob¬ 
lems  with  less  regular  geometry. 

In  the  conventional  finite  element  method,  if  a  model  with  n  degrees  of 
freedom  is  refined  by  adding  nodes,  a  new  stiffness  matrix  needs  to  be  con¬ 
structed.  In  general  the  new  stiffness  matrix  is  different  in  form  from  the 
original  stiffness  matrix.  However,  in  the  "hierarchical"  approach,  new  ele¬ 
ments  are  constructed  so  that  the  original  n  x  n  stiffness  matrix  forms  a 
submatrix  of  the  new  stiffness  matrix.  In  this  case,  therefore,  it  is  neces¬ 
sary  to  compute  only  the  stiffness  entries  corresponding  to  the  new  degrees  of 
freedom.  This  redjces  the  computing  time  in  assembly.  ( St  e  Zienkiewicz 
(ref.  6)  for  a  more  detailed  discussion  on  hierarchical  finite  elements.) 

Since  error  estimation  in  the  hierarchical  approach  has  been  established,  we 
can  take  advantage  of  the  error  analyses  in  studying  conventional  finite 
elements. 


Consider  a  finite  element  model  with  n  degrees  of  freedom  (d.o.f).  If 
the  mesh  is  to  be  refined  by  introducing  additional  nodes  (h-type  refinement), 
then  it  is  necessary  to  know  the  expected  improvement  in  error  before  a  refine¬ 
ment  step  is  undertaken.  O.C.  Zienkiewicz,  et  al.  (ref.  7)  and  Peano,  et  al . 
(ref.  8)  have  shown  that  if  the  ns-l^*1  d.o.f. 
cal  1 y ,  then  the  error  in  the  energy  norm  is: 


to  be  introduced  hierarch  *  - 


'n ,  1 


/n+l  '  Kn*l 


„0: 

,n  n  / 


(26) 


'n+1 ,n+l 
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where  fn+i  is  the  force  corresponding  to  the  n+lth  d.o.f.,  Kn+i(n+]  is  the 
stiffness  of  the  n+lth  d.o.f  Kn+]n  is  the  off-diagonal  stiffness  relating 
the  n+l^h  d.o.f.  to  the  original  n  d.o.f.  system,  and  un  is  the  array  of 
nodal  displacements  of  the  n  d.o.f.  system.  The  subscripts  n,l  of  the 
error  e  refer  to  the  n  original  d.o.f.  and  the  new  d.o.f. 

Zienkiewicz  (ref.  6)  has  used  the  above  error  relation  to  define  an  error 
indicator  in  the  form: 


(27) 


where  £  is  the  finite  element  residual. 

In  an  adaptive  refinement  strategy,  these  indicators  are  normally  calcula¬ 
ted  for  all  the  d.o.f.  corresponding  to  the  next  refinement.  The  indicators 
serve  the  purpose  of  identifying  the  region  where  refinement  is  necessary. 

Next,  the  error  corresponding  to  the  previous  iteration  wherein  the  n^h 
d.o.f.  was  added  is: 


(28) 


The  corresponding  error  indicator  is: 


(29) 


These  derivations  are  for  hierarchical  finite  elements.  However,  the 
error  with  conventional  finite  elements  will  be  of  similar  form  (ref.  6). 


The  most  general  method  of  generating  a  good  grid  is  to  have  an  equal  dis¬ 
tribution  of  some  specified  weight  function.  (See  Eiseman  (ref.  9)  for  a  com¬ 
plete  discussion  on  adaptive  grid  generation.)  Cften,  the  error  in  the  finite 
element  solution  is  used  as  the  weight  function  (ref.  10).  Therefore  the 
objective  is  to  distribute  the  error  equally  among  all  elements.  However,  the 
value  of  the  residual  £  can  be  obtained  only  after  the  equilibrium  equations 
are  solved.  Nevertheless,  one  way  of  obtaining  an  equi-di stribution  of  error 
a  priori  is  by  having  a  uniform  element  stiffness.  As  a  consequence  C,  will 
be  nearly  uniform  among  the  elements.  The  trace  minimization  procedure  deve- 


1  nnoW  Cliff 

iuj;cu  nci  c  in  yt  vuuuci  juui 


the  opti¬ 


mality  condition,  equation  (19),  is  a  constant  y 


(30) 
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Substituting  into  equation  (9),  the  element  stiffness  coefficient  is: 


which  is  a  constant.  Therefore  the  trace  minimization  procedure  produces  a 
uniform  element  stiffness. 

Finally,  observe  the  graphs  of  the  errors  in  figures  4  to  6.  As  mentioned 
earlier,  the  errors  are  more  equally  distributed  with  the  improved  mesh  while 
there  is  a  skewed  distribution  with  the  uniform  mesh.  In  tables  I  to  III,  the 
standard  deviation  of  the  within-element  rms  errors  is  lower  for  the  improved 
mesh  than  it  is  for  the  uniform  mesh.  Therefore,  in  all  cases,  the  improved 
mesh  distributes  the  error  more  uniformly  than  the  uniform  mesh.  The  rms 
errors  are  exactly  equal  (within  precision  of  calculation)  in  the  case  where 
temperatures  are  specified  at  the  boundaries.  Therefore  the  mesh  obtained  in 
this  case  is  optimal.  Similar  results,  however  are  not  obtained  in  the  case 
where  both  temperature  and  gradient  of  temperature  are  specified  because  of 
the  inability  of  FEM  to  strongly  satisfy  the  Neumann  boundary  conditions.  Nev¬ 
ertheless,  it  demonstrates  the  usefulness  of  the  trace  minimization  procedure 
in  a  priori  grid  refinement. 


CONCLUSIONS 

The  new  method  of  finite-element  grid  improvement  based  on  the  minimiza¬ 
tion  of  the  trace  of  the  stiffness  method  has  been  extended  to  the  problem  of 
heat  transfer  in  a  solid  body.  In  elasticity  problems,  this  procedure  is 
equivalent  to  minimizing  the  potential  energy  of  the  model  by  dividing  the 
strain  energy  equally  among  the  elements.  The  following  conclusions  were  made: 

1.  Nodal  positioning  obtained  by  minimizing  the  trace  of  the  stiffness 
matrix  leads  to  an  improved  mesh  over  that  obtained  by  uniform  positioning  of 
the  nodes. 

2.  Since  trace  minimization  is  an  a  priori  method,  the  mesh  may  be  refined 
without  solving  the  finite  element  problem.  This  makes  the  minimization  proce¬ 
dure  computationally  inexpensive  to  perform.  The  mesh  resulting  from  trace 
minimization  may  be  used  as  a  starting  mesh  for  other  mesh  refinement  procedures 
such  as  element  division  or  element  enhancement  (h-methods  or  p-methods). 

3.  The  method  does  not  give  satisfactory  results  for  analyses  in  which 
both  temperature  and  gradient  of  temperature  are  specified  boundary  conditions 
due  to  the  inability  of  the  FEM  to  strongly  satisfy  Neumann  boundary  conditions. 
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TABLE  I.  -  COMPARISON  BETWEEN  UNIFORM  MESH,  IMPROVED  MESH,  AND  THEORETICAL  VALUES  FOR 
TEMPERATURE  SPECIFIED  BOUNDARY  CONDITION  PROBLEM 


Radius, 

mm 

Temperature  value,  3C 

Uni  form  mesh  error 

Improved  mesh  error 

Uni  form 
mesh 

Improved 

mesh 

Theory 

At  node 

In  element 
rms 

At  node 

In  element 
rms 

20. 

25.  1487 
27.5 
31.6228 
35. 

39.7635 

'12.5 

50. 

100. 

100. 

75.0000 

100. 

75. 

65.2453 

50. 

38.0260 

n 

0.0000 

0.0000 

0.5171 

65.3550 

-0. 1096 

1.047  ) 

50.0000 

0.0000 

0.5171 

39.0247 

-0.0988 

0 .6088 

25.0000 

0.0000 

0.5171 

17.7907 

0. 

17.7  J66 
0. 

-0.0541 

0. 

0.4422 

0.2857 

0. 

0.0000 

0.5171 

Standard 

deviation  of  wi thin-element  rms  errors 
Overall  rms  error 

0.3309 

0.6801 

0.00001 

0.5210 

TABLE  II.  -  COMPARISON  BETWEEN  UNIFORM  MESH,  IMPROVED  MESH,  AND  THEORETICAL  VALUES  FOR 
TEMPERATURE/TEMPERATURE-GRADIENT  SPECIFIED  BOUNDARY  CONDITIONS 


Radius, 

mm 

Temperature  value,  °C 

Uni  form  mesh  error 

Improved  mesh  error 

Uni  form 
mesh 

Improved 

mesh 

Theory 

At  node 

In  element 
rms 

At  node 

In  element 
rms 

20. 

25.1487 

27.5 
31.6228 
35. 

39.7635 

42.5 

50. 

99.4772 

99.5700 

74.6775 

100. 

75. 

65.2453 

50. 

38.9260 

25. 

17.7360 

0. 

0.5228 

HK ! !  9 

mA  :•  9 

0.2403 

65.0133 

0.2321 

0.6786 

49.7850 

1 

0.2980 

38.8207 

0.  ’053 

0.4273 

24.8925 

0.1075 

0.3782 

17.6977 

0. 

0.0389 

0. 

0.3082 
0.24  {0 

0. 

0.0000 

0.4694 

Standard  deviation  of  wi thin-element  rms  errors  0.1920 

Overall  rms  error  0.4494 

TABLE  III.  -  COMPARISON  BETWEEN  UNIFORM  MESH,  IMPROVED  MESH,  AND  THEORETICAL 
TEMPERATURE  GRADIENT  VALUES  FOR  TEMPERATURt/TEMPERATURE  GRADIENT 
SPECIFIED  BOUNDARY  CONDITIONS 


Radius , 
mm 

Temperature  gradient,  "C/mm 

Uni  form  mesh  error 

Improved  mesh  error 

Uniform 

mesh 

Improved 

mesh 

Theory 

At  node 

In  element 
rms 

At  node 

In  element 
rms 

20. 

25.1487 

27.5 
31.6228 
35. 

39.7635 

42.5 

50. 

-4.5952 

-4.8347 

-3.84490 

-5.4568 
-4.3396 
-3.9686 
-3.4512 
-3. 1 182 
-2 . 7446 
-2.5679 
-2.1827 

0.8616 

0.6220 

-0.4947 

0.5183 

-3.4923 

0.4762 

0.6029 

-3.0578 

-0.3934 

0.4122 

-2.8164 

0.3018 

0.3594 

-2 . 43 i 7 

-0.3129 

0.3278 

-2.3597 

-2.3597 

0.2082 
-0. 1770 

C.2392 

0.1134 

-2.4317 

-0.2490 

0. 1653 

Standard  deviation  of  within-element  rms  errors  0.2086 

Overall  rms  error  0.3755 

0.1490 

0.3478 

COMPUTE  THE  STIFFNESS  MATRIX 


ERROR  IN  TEMPERATURE.  °C  TEMPERATURE, 


RADIAL  POSITION.  MM 
TEMPERATURE  B.C. 


FIGURE  3.  -  SAMPLE  PROBLEM:  THEORETICAL  AND 
NODAL  TEMPERATURE  VALUES. 


FIGURE  R.  -  EXAMPLE  1  (TEMPERATURE/TEMPORARY 
BOUNDARY  CONDITIONS):  ERRORS  IN  FEM  RESULTS 
FOR  CYLINDER  TEMPERATURE. 


TEMPERATURE/GRADIENT  B.C. 

FIGURE  5.  -  EXAMPLE  2  C  GR  AD  I  ENT/TEMPER  AT  L'RE 
BOUNDARY  CONDITIONS):  ERRORS  IN  FEM  RESULTS 
FOR  CYLINDER  TEMPERATURE. 


GRADIENT  SOLUTION  -  TEMPERATURE/GRADIENT  B.C. 

FIGURE  6.  -  EXAMPLE  2  (GRADIENT/TEHPERATURE 
BOUNDARY  CONDITIONS):  ERRORS  IN  FEM  RESULTS 
FOR  CYLINDER  TEMPERATURE  GRADIENT. 
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